#include <cmath>
#include <sstream>
#include <fstream>
#include <iostream>
#include <vector>
#include <sys/stat.h>

void exitmessage(std:: string file){
  std:: cout << "  " << std:: endl;
  std:: cout << file << " does not exsist ... " << std::endl;
  std:: cout << " check it out, for now I will STOP HERE!!! " << std::endl;
  std:: cout << "  " << std:: endl;
  exit(1);  
}

int main(){
  int iStatus = 0;
  int id0;
  int nhaloes;
  std:: vector<double> m,x,y,theta;
  std:: vector<double> rvir;

  /* structure of the info_haloes.dat */
  int n1 = 10;
  int n2 = 1;
  int n3 = 27;

  /* read mergingmasses.d */
  std:: string filea = "mergingmasses.d";
  std:: string fileb = "info_haloes.dat";
  std:: ifstream ifilea, ifileb;
  ifilea.open(filea.c_str());
  if(ifilea.is_open()){
    ifilea >> id0;
    double am, ax, ay, atheta;
    while(ifilea >> am >> ax >> ay >> atheta ){
      m.push_back(am);
      x.push_back(ax);
      y.push_back(ay);
      theta.push_back(atheta);
    }
    nhaloes = m.size();
    std:: cout << "  " << std:: endl;
    std:: cout << " ID of the main halo = " << id0 << std:: endl;
    std:: cout << " number of merging clusters " << nhaloes << std:: endl;
    std:: cout << "  " << std:: endl;
  }else{
    exitmessage(filea);
  }

  std:: vector<std:: vector<double> > a(n1,nhaloes),c(n3,nhaloes);
  std:: vector<std:: vector<int> > b(n2,nhaloes);

  ifilea.close();
  ifileb.open(fileb.c_str());
  if(ifileb.is_open()){
    for(int i=0;i<nhaloes;i++){
      for(int j=0;j<n1;j++){
	ifileb >> a[j][i];
      }
      for(int j=0;j<n2;j++){
	ifileb >> b[j][i];
      }
      for(int j=0;j<n3;j++){
	ifileb >> c[j][i];
      }
      //                Mass      and      Rvir
      // std:: cout << a[0][i] << "  " << c[21][i] << std:: endl;
    }
    ifileb.close();
  }else{
    exitmessage(fileb);
  }

  std:: cout << " I will now merge the catalogues ... " << std:: endl;

  std:: string fileg,files,sid;
  std:: ifstream ifileg, ifiles;  

  std:: vector<double> s1,s2,s3,s4,s5;
  std:: vector<double> g1,g2,g3,g4;
  s1.resize(0);
  s2.resize(0);
  s3.resize(0);
  s4.resize(0);
  s5.resize(0);
  g1.resize(0);
  g2.resize(0);
  g3.resize(0);
  g4.resize(0);
  for(int i=0;i<nhaloes;i++){
    int id = id0 + i;
    std:: ostringstream osid;
    osid << id;
    sid = osid.str();
    files = "satellites/satinfo."+ sid +".sub";
    fileg = "satellites/galaxies."+ sid +".dat";
    ifiles.open(files.c_str());
    ifileg.open(fileg.c_str());
    std:: cout << files << std:: endl;
    std:: cout << fileg << std:: endl;
    double aa,bb,cc,dd,ee;
    while(ifiles >> aa >> bb >> cc >> dd >> ee){
      s1.push_back(aa);
      // rescale with respect to the virial radius of the "central" cluster
      s2.push_back(bb*c[21][i]/c[21][0]);
      s3.push_back(cc*c[21][i]/c[21][0]);
      s4.push_back(dd*c[21][i]/c[21][0]);
      s5.push_back(ee);
    }
    while(ifileg >> aa >> bb >> cc >> dd){
      g1.push_back(aa);
      g2.push_back(bb);
      g3.push_back(cc);
      g4.push_back(dd);
    }
    ifiles.close();
    ifileg.close();
  }
  
  std:: ostringstream osid,osid2;
  std:: string sid2;
  /* file that I'm going to look for */
  std:: string filerayshoot = "rayshoot";
  std:: string filefits1 = "SHalo";
  std:: string filefits2 = "SkyLens";
  std:: string filecrit = "critical";
  std:: string fileMAPradialprof = "MAP_radial_prof";
  osid << id0;
  sid = osid.str();
  int idfinal = id0 + nhaloes - 1;
  osid2 << idfinal;
  sid2 = osid2.str();
  files = "satinfo."+ sid +".sub";
  fileg = "galaxies."+ sid +".dat";
  std:: ofstream ofileg, ofiles;  
  ofiles.open(files.c_str());
  ofileg.open(fileg.c_str());
  int ng = g1.size();
  int ns = s1.size();
  
  for(int i=0;i<ns;i++){
    ofiles << s1[i] << "  " << s2[i] << "  " << s3[i] << "   " << s4[i] 
	   << "   " << s5[i] << std:: endl; 
  }

  for(int i=0;i<ng;i++){
    ofileg << g1[i] << "  " << g2[i] << "  " << g3[i] << "   " << g4[i] 
	   << std:: endl; 
  }
  ofiles.close();
  ofileg.close();
  
  /*  bash commands */
  std:: string command;
  command = "rm -fr fitsbis satellitesbis";

  mkdir("fitsbis",S_IXOTH|S_IWRITE|S_IROTH|S_IRGRP|S_IXGRP|S_IREAD|S_IEXEC);
  mkdir("satellitesbis",S_IXOTH|S_IWRITE|S_IROTH|S_IRGRP|S_IXGRP|S_IREAD|S_IEXEC);

  command = "cd fits && ls -1 " + sid2 + filerayshoot + "*.dat 2> /dev/null > ../filelist.txt && cd ../";
  system (command.c_str());
  command = "cd fits && ls -1 " + sid2 + filefits1 + "*.fits 2> /dev/null >> ../filelist.txt && cd ../";
  system (command.c_str());
  command = "cd fits && ls -1 " + sid2 + filefits2 + "*.fits 2> /dev/null >> ../filelist.txt && cd ../";
  system (command.c_str());
  command = "cd fits && ls -1 " + sid2 + filecrit + "* 2> /dev/null >> ../filelist.txt && cd ../";
  system (command.c_str());
  command = "cd fits && ls -1 " + sid2 + fileMAPradialprof + "*.dat 2> /dev/null >> ../filelist.txt && cd ../";
  system (command.c_str());
  
  int ntype = 0;
  if(idfinal>0)  ntype = int(log10(idfinal));

  std:: ifstream fin ("filelist.txt");   
  std:: string filei,fileo,file;
  while(fin >> file){
    filei = file;
    file.erase (file.begin()+ntype); 
    fileo = sid + file;
    command = "cp fits/" + filei + " fitsbis/" + fileo; 
    system (command.c_str());    
  }

  command = "mv " + files + " satellitesbis/";
  system (command.c_str());
  command = "mv " + fileg + " satellitesbis/";
  system (command.c_str());

  command = "rm -fr filelist.txt";
  system (command.c_str());
  command = "rm -fr fits";
  system (command.c_str());
  command = "rm -fr satellites";
  system (command.c_str());
  command = "awk '{ print $18, $19} ' info_haloes.dat | tail -1 > EffAndMedThetaEmergingclusters.dat";
  system (command.c_str());
  command = "head -1 info_haloes.dat > info_haloes.dat1";
  system (command.c_str());
  command = "mv info_haloes.dat1 info_haloes.dat";
  system (command.c_str());
  command = "mv fitsbis fits";
  system (command.c_str());
  command = "mv satellitesbis satellites";
  system (command.c_str());

  std:: cout << " ... done, check it out! ..... ;-) " << std:: endl;

}
